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ABSTRACT 

Semi-analytic models of self-gravitating discs often approximate the angular momentum 
transport generated by the gravitational instability using the phenomenology of viscosity. This 
allows the employment of the standard viscous evolution equations, and gives promising re- 
sults. It is, however, still not clear when such an approximation is appropriate. 

This paper tests this approximation using high resolution 3D smoothed particle hydro- 
dynamics (SPH) simulations of self-gravitating protostellar discs with radiative transfer. The 
nature of angular momentum transport associated with the gravitational instability is char- 
acterised as a function of both the stellar mass and the disc-to-star mass ratio. The effective 
viscosity is calculated from the Reynolds and gravitational stresses in the disc. This is then 
compared to what would be expected if the effective viscosity were determined by assuming 
local thermodynamic equilibrium or, equivalently, that the local dissipation rate matches the 
local cooling rate. 

In general, all the discs considered here settle into a self -regulated state where the heating 
generated by the gravitational instability is modulated by the local radiative cooling. It is found 
that low-mass discs can indeed be represented by a local a-parametrisation, provided that the 
disc aspect ratio is small (H/ R ^ 0.1) which is generally the case when the disc-to-star mass 
ratio q < 0.5. However, this result does not extend to discs with masses approaching that of 
the central object. These are subject to transient burst events and global wave transport, and 
the effective viscosity is not well modelled by assuming local thermodynamic equilibrium. 
In spite of these effects, it is shown that massive (compact) discs can remain stable and not 
fragment, evolving rapidly to reduce their disc-to-star mass ratios through stellar accretion 
and radial spreading. 
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1 INTRODUCTION 

Accretion discs play an important role in many astrophysical sit- 
uations, from protostellar systems through to discs around super- 
massive black holes in active galactic nuclei (AGN). What is still 
very uncertain is the process through which angular momentum is 
transported outwards in such discs. It is clear, from observations 
of accretion rates, that classical hydrodynamical viscosity is in- 
sufficient to play this role. The typical solution is to assume an 
ad hoc parametrisation of the viscosity, whose origin is not well 
understood. The archetype is the a-parametrisation introduced by 
IShakura & Sunvaevl i fl973h in which the viscosity v is assumed to 
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depend on the disc sound speed, c s , and thickness, H, through 
v — ac s H, where a « 1. 

This allows a number of different physical mechanisms to be 
considered as the origin of this viscosity. The most frequently in- 
voked is turbulent viscosity, shifting the problem to the origin of 
the turbulence. If t he disc is sufficiently ionised, the magnetoro- 
tational instability jBalbus & Hawlevlll99ll : fealbus & Papal oizoul 
ll999l : |Papal oizou & Nelson 120031) can result in turbulence that can 
provide the required viscosity. However, if the disc is very weakly 
ionised (as in the case of most protostellar discs at early times), 
another source must be sought. During the earliest stages of star 
formation, when disc masses are likely to be high relative to the 
mass of the central protostar, gravitational instabilities (GI) may 
provide the answer (Lin & Pringle 1987; Laughlin & Bodenheimer 
Il994l) . 
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The susceptibility of an infinitesimally thin disc to gravita- 
tional instabili ty can be measured using the Toomre Q parameter 
llToomrell 19641) : 



Q 



ttGE' 



(1) 



where c s is the local sound speed, E is the disc surface density and 
k is the epicyclic frequency (equal to the angular frequency fl in 
Keplerian discs). Discs are gravitationally unstable to axisymmetric 
ring perturbations if Q < 1, while simulations have shown that for 
Q < 1.5 — 1.7 discs are unstable t o the growth of nonaxisymmetric 
perturbations (Durisen et al. 2007). Gravitational instabilities will 
generally lead to a self-regulating, quasi-steady state in such discs 
(Paczvnski 1978). Discs that are cool enough to become unstable 
will be heated by the gravitational instability through shocks, in- 
creasing their Q until they reach stability. Discs that are initially 
too hot for the instability to set in will undergo radiative cool- 
ing towards instability. These competing processes control the disc 
thermodynamics such that the value of Q is kept close to, but just 
above , the instability boundary and is referred to as marginal sta- 
bility ( |Paczvnskilll978l:lBertin&Lodatolll999h . 

However, to put forward turbulence generated by gravitational 
instability as the source of the unknown "viscosity", the nature of 
the angular momentum transport generated in this manner must 
be investigated. In particular, can the Q-parametrisation be used 
to evaluate the viscosity generated by the gravitational instabil- 
ity? If this approxi mation is to be used, then th e transport needs 
to be local in origin. iBalbus & Papaloizoul ( fl99gh have shown that 
the energy flux generated by GIs contains terms that are inherently 
non-local (associated with global wave transport), indicating that 
the phenomenology of viscosity will never exactly reproduce the 
tra nsport induced by grav itational instabilities. However, as shown 
by lLodato & Ricel d2004l) . the a-approximation may be sufficient 
to explain disc behaviour in systems where global wave transport 
is negligible. Therefore, the problem can be addressed by consider- 
ing some key questions: 

Is angular momentum transport local? Can an effective viscous a 
be estimated from the assumption of local thermodynamic equi- 
librium ? Do realistic, self-gravitating protostellar discs settle into 
marginally-stable, quasi-steady states ? 

Previous work on the locality of this angular momen- 
tum transport has relied heavily on numerical simulations. 
lLaughlin & Rozvczk j Jl996l) used 2D grid based simulations to in- 
dicate that the value of a must vary with orbital radius (to produce 
th e expected density evolution). In three dimensions, the early work 
of Laughlin & Bodenheimerj j 19941) using smoothed particle hydro- 
dynamics (SPH) simulations of massive, isothermal discs showed 
that simple a models do indeed reproduce the correct density evo- 
lution. However, the strength of the gravitatio nal instability is in- 
herently linked to the disc thermodynamics jPickett et alj l200(j| ; 
iNelson et al]|2000l) . Any physically realistic study of angular mo- 
mentum J^msrjortjw sejfjjjrayky^ include radiative 
effects jPickett et al.l 120031; iMeiia et al.ll2005t). F ollowing the ap- 
proach of iGammiel fcOOll) . lLodato & Ricel fe004l) used SPH simu- 
lations with an adiabatic equation of state, but with a cooling time 
of the following form 



tcooi^ = f3 — constant. 



(2) 



With the above cooling, the local approximation would suggest that 
jGammiejl200lh . 



(3) 



lLodato & Ricel ^2W^) show that this approximation is valid, 
and that transport is local, in discs with mass ratios q = Md/M 4 
less than 0.25 (and aspect ratios H/R ^ 0.1), where the self- 
regulation controlled by Q ensures a quasi-steady state. Fu rther in- 
vestigation of more massive discs jLodato & Ricel 120051) showed 
that despite the evolution being clearly non-steady (with recurrent 
episodes of variable accretion) th ere was no significan t evidence 
for global wave energy transport. ICossins et al.1 d2009h have also 
carried out a detailed analysis of the gravitational instability under 
this cooling time approximation, investigating discs with q < 0.1 
and char acterising the resultant spiral structure. They demonstrated 
(see also lBalbus & Papaloizoul \l 9991) that global transport occurs 
whenever spiral waves dissipate far fr om their corot a tion r adius. 
For the low-mass ratios considered in ICossins et alj d2009l) . this 
does not happen and the resulting transport is therefore local and 
quasi-steady to a high degree, showing that the viscous approxima- 
tion works for discs with parametrised radiative cooli ng, although 
it may depend on the form of the cooling function (Mejia et al. 
| 2005 |: iDurisen et alj 120071). Recent semi-ana lytic works faarkel 
l20091 : lRice & Armitagell2009l ; lRice etZll20ld) has, however, used 
this approximation to study the formation and evolution of massive 
protostellar discs. 

This paper builds on these earlier results using global 3D SPH 
simulations of protostellar discs over a range of stellar masses and 
disc-to-star mass ratios. What makes this different to most earlier 
work is that the SPH code, in this c ase, uses a hybrid method of 
radiative transfer (Forgan et al. 2009), which models the effects of 
frequency-averaged radiative transfer without significant runtime 
losses. By adding radiative transfer, these simulations are in the 
best position to accurately model gravitational instabilities in re- 
alistic protostellar discs. The analysis will focus on the key ques- 
tions defined earlier, in effect to characterise the efficacy of the q- 
parametrisation in self-gravitating protostellar discs. Section[2]will 
outline the key physics involved in this work; section [3] will focus 
on the numerical techniques used to produce the simulations; sec- 
tion [4] will outline and discuss the results of the simulations, and 
section[5]will summarise the work. 



2 ANGULAR MOMENTUM TRANSPORT AND THE 
a -PAR AMETRIS ATION 

If a thin disc approximation is adopted, an accretion disc's equa- 
tions of motion can be cast in terms of vertically averaged proper- 
ties. Therefore, the equation of continuity (using cylindrical coor- 
dinates) becomes 



as 1 d . ^ . „ 



(4) 



where E is the surface density that depends on position, r, and time, 
t, and v-c is the radial velocity of the disc material. Conservation of 
angular momentum gives 



| (Er 2 n) + i£ (Er 3 fK) =l± (r 2 T r ,) 
at v ' r or v ' r or v ' 



(5) 



where T It f, is the (vertically averaged) viscous stress tensor compo- 
nent. The calculation of T T( p is the crux of the problem, and the most 
important facet of accretion disc theory in general. As has already 
been stated, typical hydrodynamical v iscosity is insufficient. To 
characterise T x <j>, the Q-parametrisation ( Sha kura & Sunvaevl 19731) 
can be used: 



) 7(7 - l)i C oolfi 
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dlnr 



(6) 
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or equivalently, in terms of the kinematic viscosity v. 



3 METHOD 



v — ctc B H, 



(7) 



where H — Cs/fi is the scale height of the disc. If the disc is in 
thermal equilibrium, an expression can be found for a by equating 
the rate at which viscosity dissipates energy in the disc with the 
rate at which this energy is lost through radiative cooling. Viscous 
dissipation occurs according to 



Q 



(8) 



where this describes the dissipation rate per unit surface. The radia- 
tive cooling can be parametrised in terms of the local cooling time, 
*cooi> giving 



cr 



V 
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(9) 



where U is the internal energy per unit surface, and 7 is the ratio 
of specific heats. Equating Q + and Q~ and rearrangin g gives the 
following expression for a jPringldl 1 98 ll : lGammiell200 II) 



(10) 



dlnr J 7(7 — l)t coo if2 

Note that equation ( |10t requires that local heating and cooling be 
in balance: in practice, this balance must be true over some char- 
acteristic timescale, where we should instead equate time averaged 
quantities, i.e. < Q + >~< Q~ >. If the disc is self-gravitating, 
the component of the viscous st ress tensor associated wi th the grav- 
itational instability is given by l ILvnden-Bell & Kalnaisll 19721) 



± v<t> — 



9r94- 



dz, 



(11) 



where g T and are the components of the gravitational accelera- 
tion in cylindrical coordinates. The full viscous stress tensor also 
includes the 'Reynolds' stresses (i.e., stresses produced by velocity 
and density perturbations as a result of gravito-hydrodynamics) 



(12) 



where 8v T and 61)4, are (vertically averaged) fluctuations from the 
mean fluid velocity (again in cylindrical coordinates). The total vis- 
cous stress in the disc is therefore the sum of these two tensor 
components. Using T vt j, = T^ yn + T^ av together with equation 
(|6} provides a means for calculating an effective a associated with 
gravitational instabilities 



Qtotal 



dlnQ, 
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If the angular momentum transport is local, the stress tensor, 
and consequently atotah depend only on local conditions in the 
disc and equation dlOt would also be valid. Gravitational stresses 
may, however, be exerted as a result of global features in the po- 
tential at large se parations (such as spiral den sity waves). In fact, it 
has been shown (Balbus & Papaloizou 1999) that the energy trans- 
port associated with gravitational instabilities contains global terms 
and, if such terms are significant, a local prescription for angular 
momentum transport in self-gravitating discs may be a very poor 
approximation. A prime goal of this work is to compare atotai, 
computed as above from the Reynolds and gravitational stresses in 
the disc, with a C ooi, computed by assuming that the disc is in local 
thermodynamic equilibrium. 



3.1 SPH and the Hybrid Radiative Transfer Approximation 

Smoothed Particle Hydr o dynamics (SPH ) jLucvl 1 19771 : 
iGingold & MonagharJ 1 19771 : iMonagharJ Il992h is a Lagrangian 
formalism that represents a fluid by a distribution of particles. Each 
particle is assigned a mass, position, internal energy and velocity. 
State variables such as density a nd pressure are then ca lculated 
by interpolation (see reviews by |Monaghan| [l992. 2005). In the 
simulations presented here, the gas is modelled using 500,000 SPH 
particles while the star is represented by a point mass particle onto 
which gas particles can acc rete, if they are sufficiently close and 
are bound JB ate et al .11 19951) . 

The SPH cod e used in this w ork is based on the SPH 
code developed by iBate et al.l j 19951) which uses individual par- 
ticle timesteps, and individually variable smoothing lengths, hi, 
such that the number of nearest neighbours for each particle is 
50 ± 20. The code uses a hybrid method of approximate ra- 
diative transfer dForgan et al.l I2009T) . which is built on two pre- 
existing radiati ve algorithms: the poly tropic cooling approxima- 
tion devised bvlStamatellos et alJd2007l). and flux- lim ited diffusion 
(e.g. IWhiteh ouse & B atel2004l : lMaver et al.l2007l, see lForgan et al.l 
2009 for details). This union allows the effects of both global cool- 
ing and radiative transport to be modelled, without imposing extra 
boundary conditions. 

The opacity and temperature of the gas is calculated using a 
non-trivial equation of state. This accounts for the effects of H2 dis- 
sociation, H° ionisation, He and He + ionisation, ice evaporation, 
dust sublimation, molecular absorpti on, bound-free and free-free 
transitions and electron scattering jBeH & LirJ[l994l ; I Bolev et al.l 
l2007l : IStamatellos et alj|2007l) . Heating of the disc is also achieved 
by P dV work and shock heating . 



3.2 Initial Disc Conditions 

The gas discs used in this work were initialised with 500, 000 
SPH particles located between ri n = 10 au and r out = 50 
au, distributed such that the initial surface density profile was 
E oc r~ 3 ^ 2 and with an initial sound speed profile of c s oc 
r _1 / 4 \Y e are primarily interested in considering quasi-steady self- 
gravitating systems, rather than systems that could fragment to 
form bound companions. These initial conditions (in particular the 
small disc radii) were therefore motivated by recent work suggest- 
ing that massive discs will fragment at radii beyond ~ 60 — 70 au 
( Rafik ovll2003: IStamatellos et al.ll2007l: IStamatellos & Whitworthl 
l2008l : IClarkell2009l : lRice & Armitagell2009l) . This result is consis- 
tent with observ ations that massive di scs tend to have outer radii 
less than 100 au (Rodriguez et al. 2005) and with observations sug- 
gesting t he presence of a pro toplanet at ~ 65 au in the disc around 
HL Tau dGreaves et al.ll2008h . A summary of the disc parameters 
investigated can be found in TableQ] The simulations were selected 
to evaluate the a-approximation's ability to function under 



(i) increasing disc-to-star mass ratio, q, and 

(ii) increasing stellar mass, M* 

As we are interested in q, which will evolve as the star accretes 
from the disc, we should be rigorous and also define qi n i t as the 
value of q at the start of the simulation. 
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Table 1. Summary of the disc parameters investigated in this work. 



Simulation 


Af, (M ) 


Sfimt = M d /M* 


M d (M Q ) 


1 


1.0 


0.25 


0.25 


2 


1.0 


0.5 


0.5 


3 


1.0 


1.0 


1.0 


4 


1.0 


1.5 


1.5 


5 


0.5 


0.25 


0.125 


6 


2.0 


0.25 


0.5 


7 


5.0 


0.25 


1.25 


8 


0.5 


1.0 


0.5 


9 


2.0 


1.0 


2.0 



3.3 Resolution 

There are several resolution requirements that mus t be discussed at 
this p oint. The first is the standard Jeans criterion dBate & Burkerj 
1 19971) . As some of the discs used in this work are very massive 
compared to the mass of the parent star, the possibility of fragmen- 
tation exists. To ensure that potential fragmentation is resolved, the 
minimum Jeans mass resolvable (one neighbour group of SPH par- 
ticles, around 50 in the case of the code used) must be sufficiently 
small: 

N neigh 



2N ncigh mi = 2M to 



N u . 



(14) 



The minimum Jeans mass resolvable ranges between 50M® for the 
most massive disc and 4A#e for the least massive. As it is expected 
that fragment masses will be typically severa l orders of magnitude 
higher than these values jKratter et al.ll201oT) . this establishes that 
the simulations would comfortably resolve disc fragmentation if it 
were to occur. 

Perhaps more important are the resolution issues raised by ar- 
tificial viscosity. While required by the SPH code used, this artifi- 
cial viscosity must be quantified so that we know where in the disc 
the artificial viscosity is likely to be lower than the effective viscos- 
ity generated by the gravitational instabiliti es. The linear term for 
the artificial viscosity can be expressed as jArtvmowicz & Lubowl 



the artificial viscosity can be expressed as ( 
ll994lMurravlll99dlLodato & Pricefe OlO) 



10 



OtSPHC s h, 



(15) 



where c s is the local sound speed, h is the local SPH smoothing 
length, and qsph is the linear viscosity coefficient used by the 
SPH code (taken to be 0.1). We can define an effective a param- 
eter associated with t he artificial viscosity by using equation 
dLodato & Rlcell2004l) 



^art 



rtC s //, 



and hence combinin ; 
dArtymowicz & Lubowl 1199 
12010) 



(16) 

equations l |15t an d dl6|> gives 
iMurravl Il99l : iLodato & Pricel 



1 h 

To aspH 



H 



(17) 



This shows that where the vertical structure is not well resolved 
(i.e., i is large), artificial viscosity will dominate. In the simula- 
tions presented here, this is likely to be the case inside ~ 10 au, so 
any data inside this region can not be used. We therefore did not ini- 
tially populate the region inside 10 au and although particles will 
move inside 10 au during the course of the simulations, we only 
consider results outside this radius. 



4 RESULTS AND DISCUSSION 

All of the simulations presented here were evolved for 27 outer ro- 
tation periods (ORPs)q This ensures that all our simulations have 
sufficient time to settle into quasi-steady states. In fact, the dura- 
tion of these simulations (~ 10 4 years) is roughly 10% of the main 
infall phase, during which we expect protostellar discs to be self- 
gravitating, and therefore we capture a significant fraction of the 
self-gravitating history of such discs. 

We consider two free parameters, the disc mass Md, and the 
disc mass ratio, q = Ma/M». Both q and the local sound speed 
determine whether a disc is self-gravitating or not. The sound speed 
is determined by the local radiative physics, in particular the optical 
depth to the midplane. The optical depth is a function of the disc 
surface density, E, which in turn is related to the disc mass, Md. It 
can then be seen that the values of both q and Ma will affect the 
disc's evolution under self-gravity. 

Secondly, there is the issue of how to calculate a C ool- The ra- 
diative transfer algorithm allows the calculation of t coo \ for each 
SPH particle, and therefore each particle has its own a C ooi- How- 
ever, equation d 1 Ot shows that particles with short cooling times 
(e.g., those at higher elevation from the midplane) can skew at- 
tempts to create azimuthalry averaged radial profiles. Therefore, 
when comparing a coo i with atotal, two quantities are considered: 
Q coo i, using the midplane values of t coo i, fl and 7, and a coo i cal- 
culated using vertically averaged values of F CO oi. and 7. We cal- 
culate t C ooi by first averaging the specific internal energy u and its 
rate of change ii separately, giving 

u 

^cool = — • (18) 

u 

This distinction between midplane and vertically averaged values 
is important. Using the midplane values of t coo \ allows us t o deter- 
mine t he v alidity of recent ID semi- analytic models, such as Clarke 
(2009) and iRice & Arm itage J2009h . that calculate transport prop- 
erties based on the midplane temperature. The vertically averaged 
quantities, however, give a more accurate estimate of the rate at 
which the disc loses energy and allows us to establish if local heat- 
ing and cooling is in balance. This will then determine if the local 
a-approximation is still appropriate, even if using midplane values 
is not. 



4.1 The Influence of Disc Mass 

To study the effect of increasing disc mass on angular momentum 
transport, Simulations 1, 2, 3 & 4, which share the same stellar 
mass (M* = IMq) are analysed together. These discs have initial 
masses of 0.25, 0.5, 1.0 and 1.5 Mq respectively. 



4.1.1 General Evolution 

Despite all four simulations beginning with a wide range of disc 
masses, their surface density profiles do not differ greatly between 
r ~ 20 — 60 au, as can be seen in Figure [2] The higher mass 
discs (ginit = 1 & (fruit = 1.5) are in general much denser be- 
tween r ~ 10 — 20 au, indicating mass bui ld-up in the inner re- 
gions as suggested and seen by ot her authors ( Armitag e et alj200ll : 
IZhu et alj|2009l : IRice et alJ[2oToh . The lower-mass discs ((/mit — 
0.25 & ginit = 0.5) undergo a period of quiescent settling lasting 



1 Outer rotation periods are defined as the rotation period at the initial outer 
radius of the disc, r ou t = 50 au, with 1 ORP equal to 354 years. 
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Figure 1. Images showing the surface density structure of Simulations 1 (top left), 2 (top right), 3 (bottom left) & 4 (bottom right) after 27 ORPs. The stellar 
mass in each case is 1 Mq, and the initial disc masses of 0.25 Mq, 0.5 Mq, 1 Mq and 1.5 Mq respectively. The axis ranges are shown in each figure and it 
is clear that the more massive discs exhibit higher amplitude spiral structures, in particular the m = 2 mode. 



approximately 2000 years, adjusting themselves by accretion onto 
the central star, spreading in radius (see Figure [TJ and by cooling 
towards marginal in stability, ultimately se ttling into quasi-steady, 
self-regulated states dLodato & Ricd l2004), 



The higher mass discs (ginit = 1 & (ftnit = 1.5) undergo sev- 
eral transient burst events, marked by persistently strong m = 2 
spiral activity (see FigureQJ. They also adjust their q more rapidly 
compared to the two lower mass discs, with reductions between 20- 
30% over approximately 10 ORPs. This is due to significant accre- 
tion, with the central star accreting a total of O.23A/0 for q^t = 1 
and 0.38A/ W for <? illit = 1.5, and is con sistent with the sugges- 
tion jRice & Arm itage 20oH lClarkd f2009) that the mass accretion 
rate has a very strong dependence on surface density or, equiva- 
lently, disc mass. The discs with q^^t > 0.5 also spread to a much 
larger radius than the gi n j t < 0.5 discs (which is clear in Figure 
[]}, with significant fractions of mass outside 60 au. All the discs in 
Figure [2] are stable against f ragmentation, with j3 = t coo \Sl » 3 
(Qcooi < 0.06) at all radii jGammidl200ll : iRice et alj|2003l) . The 
values of ft as a function of opacity regime are als o in good agree- 
ment with those predicted by Cossin set al.1 d2010h . 



Considering the azimuthal Fourier modes of the higher mass 
discs (Figure [3} confir ms previous results regarding mode strength 
and d isc mass ratio dLodato & Ricd 12004 , 120051 ; ICossins et al.l 
2009). The lower mass ratio discs have power distributed over a 
range of modes (up to m ~ 8) with the m = 2 mode (and its 
harmonics) becoming dominant as q increases, indicating the pos- 
sibility of global transport in the discs. The gi n i t = 1 disc appears 
to have a larger m — 2 amplitude than the qi n i t = 1.5 disc. The 
precise reason for this is difficult to ascertain based on the available 
evidence, but it may be due in part to a) the more rapid evolution of 
q in the latter case (Figure [T] lower right panel), and/or b) a more 
efficient cascade of power into the harmonics m = 4, 6, 8 reducing 
the amplitude at m = 2. 



4.1.2 The a Approximation 

What we really want to establish is whether or not these discs obey 
the local viscous approximation. If they do, then the effective a 
parameter for these discs can be approximated using equation J 1 0b _ 
Figure|4]shows the azimuthally averaged, radial a profiles for the 4 
simulations in which M* = IMq . The radial profiles in each case 
are also time averaged over the final 13 ORPs. In each panel, the 
solid line is atotai computed using Equation l !13b , while the dashed 
line the midplane a coo i and the dotted line a vertically averaged 

In the low-mass case ((fruit = 0.25), it can be seen (Figure[4]l 
that Qcooi calculated from both the midplane cooling time (dashed 
line) and the vertically averaged cooling time (dotted line) approx- 
imates well a to tai, computed directly from the Reynolds and grav- 
itational stresses. That atotai increases with radius beyond 15 — 20 
au is also consistent with numerical and semi-analytic calculations 
that use the local approximation for calculating the effective grav- 
itatio nal viscosity dZhu et alJ[2009l : lRice & Armitagd] 20091 : Iciarkel 
2009). The same is true for <fr n it = 0.5, but it can be seen that this 
approximation fails for the higher-mass discs in Simulations 3 and 
4, with the profile for atotai being quite different to that for the mid- 
plane a coo i. The vertically averaged q coo i is a slightly better match 
to atotai, however, the radial profiles are quite different with atotai 
being flatter than a C ooi. This shows that, for the higher-mass discs, 
the local torque - in a time-averaged sense - is different to what 
would be expected if the effective viscous dissipation rate matched 
the local cooling rate and s uggests the presence of non-local en- 
ergy transport JCossins et alj|2009l) . That atotai exceeds the verti- 
cally averaged a coo i at small radii (r < 40 au), and is less than the 
vertically averaged a CO oi at larger radii (r > 40 au) suggests that 
energy is being transported, via global wave modes, from the inner 
to the outer disc. 

Note that both of the high-mass simulations have disc aspect 
ratios above 0.1 a cross their entire disc radius, suggested to be 
a critical val ue by Lodato & Ric^ d2004l) for deviations from lo- 
cal transport. IKratter et al.1 J2008) have suggested that there should 
be two self-gravitating a parametrizations, one for when high- 
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Figure 2. Azimuthally averaged radial profiles from the M* = IMq simulations (Simulation 1 (solid line), Simulation 2 (dotted lines), Simulation 3 (dashed 
lines) and Simulation 4 (dot-dashed lines)) after 27 ORPs. The figures show the time average of each variable (taken from the last 13 ORPs, to give the discs 
time to settle into quasi-steady states). The top left panel shows the surface density profile, the top right shows the aspect ratio, the bottom left shows the 
midplane temperature, and the right hand panel shows the disc-to-star mass ratio, q, as a function of time. Artificial viscosity dominates inside 10 au, so data 
from inside this region is ignored. 
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Figure 3. Azimuthal mode amplitudes for the M* = IMq simulations (Simulation 1 top left, Simulation 2 top right, Simulation 3 bottom left, Simulation 4 
bottom right). The figures show the time average of the modes (taken from the last 13 ORPs). These figures illustrate how the m = 2 mode becomes more 
dominant as the disc-to-star mass ratio, q, increases, indicating the presence of large-scale, global spiral density waves. 
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Figure 4. Azimuthally averaged a parameter, time-averaged over the last 13 ORPs of the simulations (Simulation 1 top left, Simulation 2 top right, Simulation 
3 bottom left and Simulation 4 bottom right). The solid line indicates the a calculated from Reynolds and gravitational stresses, the dashed line indicates 
Q cool calculated using the midplane cooling time, while the dotted line indicates a coo i calculated from the vertically averaged cooling time. For illustrative 
purposes, we also show the stress tensor component due to gravitational instability « gr av, indicated by the dot-dashed line, and the the stress tensor component 
due to the artificial viscosity a ar t . 



m modes dominate and another for when low-m modes domi- 
nate. Our results would suggest that there is some merit in this 
suggestion with the local approximation being appropriate when 
9init < 0.5, changing to an approximately radially independent a 
when ginit > 0.5. Fixing the value of a in the latter case appears 
difficult although our results may suggest that the value derived 
from the local approximation at r ~ 40 au may be suitable. 

The increase of a with decreasing radius inside 20 au is a re- 
sult of the numerical visocity cwt (the triple-dot dashed lines in 
Figure dominating in these inner regions, illustrating why we 
do not consider the region inside 10 au. The dash-dot lines in Fig- 
ure [4] show the effective gravitational a computed using only the 
gravitational stresses (i.e., a gra v = (dlnr2/dlnr) _1 T^ av /Ecf). 
This illustrates that in the inner disc, due to the dominance of the 
numerical viscosity (triple-dot dashed lines), the Reynolds stresses 
dominate over the gravitational stresses. If we were able to reduce 
th e numerica l visco sity significantly we would ex pect, as suggested 
bv lZhu et alj j2009h and lRice & Armitagd d2009h . that the effective 
gravitational a in the q < 0.5 simulations would continue decreas- 
ing to very small values in the inner disc, potentially leading to a 
pile-up of mass and periodic FU Orionis-like outbursts if the tem- 
perature in the these inner regions becomes high enough for MRI 
to operate jArmitage et al.ll200ll ; lz"hu et al.ll2009h - 

4.1.3 Are the discs quasi-steady? 

Although the mismatch between the a t otai profiles and the a coo i 
profiles in the higher-mass simulations (see Figure [4]l suggests the 
presence of non-local transport, it does not tell us whether these 
simulations reach quasi-steady states or not. To identify how quasi- 
steady the discs are, the discs' temperature profiles and Toomre in- 



stability profiles are averaged over the final 13 ORPs. The standard 
deviation about this mean is then measured, and the normalised 
quantities o t JT and uq/Q are calculated for each radius (Figure 
O. This shows the deviation of the disc from quasi-steady, thermal 
equilibrium (through <jt /T) and its deviation from a marginally- 
stable, self-regulated state (through <jq/Q). 

Simulation 1 ((fruit = 0.25, solid line in Figure [5} shows 
the lowest temperature deviation, maintaining thermal balance to 
within around 5% except in the outer regions, where this is mainly 
due to the reduced value of T. A deviation of 1 K from a mean 
of 10 K will be more significant than from a mean of 100 K. This 
is also true for g init = 0.5 (dotted line in Figure[5](, although the 
amplitude increases further at larger radii. The lower-mass simula- 
tions (gtnit < 0.5) are therefore not only local, but also settle into 
long-lived, quasi-steady states. 

The temperature profiles for the high-mass (ffrnlt > 0.5) discs 
(dashed and dash-dot lines in Figure O show significant variation 
(varying by as much as 60% in the worst case), illustrating that 
these discs not only have non-local transport, but also do not at- 
tain well-defined, long-lived quasi-steady states. This implies that 
in these discs - at any given location - there will be periods when the 
dissipation rate exceeds the local cooling rate (causing the tempera- 
ture to rise) followed by a period when the cooling rate dominates. 
This is presumably inherently linked to the global nature of the 
energy transport in these simulations. Energy is being transported 
non-locally, and is hence not being generated and dissipated at the 
same location, and therefore it is not possible for the local heating 
and cooling rates to balance at all locations in the disc. 

Figure [5] also shows deviations from uniform Q, with again 
the lower-mass discs showing the lowest deviation in the inner 50 
au, averaging around 10%. Simulations 3 & 4 (ginit > 0.5) again 
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Figure 5. Variation in the mean temperature profile (left panel) and the mean Toomre instability profile (right panel) for the M t = IMq simulations 
(Simulation 1 (solid lines), Simulation 2 (dotted lines), Simulation 3 (dashed lines) and Simulation 4 (dot-dashed lines)), averaged over the last 13 ORPs. 



vary much more significantly, peaking at around 40%. These results 
show that for qinit > 0.5 a disc is unable to settle into a long-lived, 
marginally-stable, self-gravitating state. 



4.1.4 Is the transport non-local? 

Although the above suggests that there is non-local transport in the 
higher-mass discs, we have not yet convincingly shown that this is 
indeed the case. One way to do this is to compare the pattern speed 
of the dominant spiral mode, Q p , with the angular speed of th e 
disc material itself, fi. As shown bv lBalbus & Papaloizoul i 19991) . 
transport through gravitational instability can only be described in 
viscous terms when fl p — Q. When fl p ^ Q, the non-local trans- 
port terms become more significant. Waves producing non-local 
transport therefo re have a pattern speed that deviates significantl y 
from corotation (Balbus & Papaloizou 199^: ICossins et al", I l2009h . 
Equivalently, the non-local transport f raction £ must deviate signif- 
icantly from zero dCossins et al. 2009]), 1 



, where 



(19) 



j can be calculate d from the dispersion relation fo r finite 
thickness Keplerian discs (Bertin 2000; Co ssins et all 2009) 



m (Q p — fl) = c B k 



27rGE|fe| 2 
l + \k\H + 



(20) 



The factor of 1 + \k\H is required as the disc thickness dilutes 
the vertical gravitational potential. In order to determine the dom- 
inant modes, the radial and azimuthal wavenumbers (k, m) are 
spectrally averaged for each radius (i.e., the average is weighted 
by the squared amplitude in each mode), and hence Q, p is calcu- 
lated for each radius, which allows the calculation of £(r), shown 
in Figure [6] (where we have averaged £ over the last 13 ORPs). As 
can be seen, £ increases with increasing disc mass, exceeding 1 for 
Omit J? 1, illustrating that non-local transport becomes important 
as the disc-to-star mass ratio exceeds 0.5. The most massive disc 
(?init = 1-5) undergoes rapid evolution to adjust its q towards 0.85 
with a flat surface density profile, ensuring that £ is also flat out to 
larger radii (exceeding the ginit = 1 disc outside 40 au). The peak 
values of £ at around 20 — 30 AU are consistent with the peak de- 
viations of atotaj from Qcooii lending weight to the conclusion that 
non-local effects transport energy from the inner disc to the outer 
disc. 



q=C 25 - 
q = 0-5 ■ 



r (AU) 

Figure 6. The non-local transport fraction, £, for the Mr = IMq simula- 
tions (Simulation 1 (solid lines), Simulation 2 (dotted lines), Simulation 3 
(dashed lines) and Simulation 4 (dot-dashed lines)), averaged over the last 
13 ORPs. 



4.2 The Influence of Stellar Mass 

To disentangle the influences of disc mass and disc-to-star mass ra- 
tio, two sets of simulations are to be analysed together. The first 
set of discs have q- ln i t = 0.25 (Simulations 1, 5, 6 & 7), but 
have different stellar mass. The previous section showed that the 
Q-approximation holds well for Simulation 1. If disc-to-star mass 
ratio is the key property that determines the nature of angular mo- 
mentum transport (and not the local sound speed), then the a- 
approximation should be equally effective for all simulations in this 
first set. 

The second set will analyse the discs with qx n i% = 1 (Sim- 
ulations 3, 8 & 9). If q is key to the nature of angular momentum 
transport, then it should be expected that non-local transport should 
be exhibited by all the discs in the second set. 



4.2.1 The qinit = 0.25 discs 

4.2.1.1 General Evolution As with the previous set of simu- 
lations, the discs undergo an initial settling phase, and become 
marginally-stable after a period of cooling. The low initial value 
of q is relatively unchanged in all simulations, with the most mas- 
sive disc changing mass by less than 20% (see Figure [8] bottom 
right panel). All four simulations share similar aspect ratios - this 
follows from the result that the aspe ct ratio H/r is proportional 
to q during marginal instability (c.f., lLodatoll2008h . For this to be 
possible, the surface density profiles must therefore increase with 
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Figure 8. Azimuthally averaged radial profiles from the <j init = 0.25 simulations (Simulation 5 (solid line), Simulation 1 (dotted lines), Simulation 6 (dashed 
lines) and Simulation 7 (dot-dashed lines)) after 27 ORPs. The figures show the time average of each variable (taken from the last 13 ORPs). The top left 
panel shows the surface density profile, the top right shows the aspect ratio, the bottom left shows the midplane temperature, and the right hand panel shows 
the disc-to-star mass ratio, q, as a function of time. Artificial viscosity dominates inside 10 au, so data from inside this region is ignored. 



disc mass, as can be seen in the top left panel. However, the radial 
dependence of the surface density is roughly the same for all discs. 
This in turn ensures the more massive discs are hotter (bottom left 
panel), with similar radial temperature profiles for all four simula- 
tions. 



4.2.1.2 The a Approximation Repeating a similar analysis of 
a as was done for the M* = IMq discs, it can be seen (Fig- 
ure^ that the Q-approximation holds with increasing stellar mass, 
confirming that the key parameter is the disc-to-star mass ratio, q, 
which is held constant here. A local approximation therefore ap- 
pears to be suitable for systems in which qinit < 0.5. Simulation 
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Figure 9. The a parameter for the qi n i t = 0.25 simulations (Simulation 5 top left, Simulation 1 top right, Simulation 6 bottom left & Simulation 7 bottom 
right), averaged over the last 13 ORPs of the simulations. The black line indicates the a calculated from the Reynolds and gravitational stresses (atotalX me 
dashed line indicates o coo i calculated using the midplane cooling time at that radius, and the dotted line indicates the a coo i calculated from the vertically 
averaged cooling time. 
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Figure 10. Variation in the mean temperature profile (left) and the mean Toomre instability profile (right) for the q = 0.25 simulations (Simulation 5 (solid 
line), Simulation 1 (dotted lines), Simulation 6 (dashed lines) and Simulation 7 (dot-dashed lines)), averaged over the last 13 ORPs. 



7 in which M* = 5M© suggests that there may be some depen- 
dence on the stellar mass as the calculated atotai is somewhat lower 
than the expected q coo i inside 30 au. The aspect ratio of this disc 
is, however, quite flat with H/r > 0.1 for a much wider radial 
range than in the other simulations. The region where the aspect 
ratio exceeds 0.1 corresponds with the region where atotai devi- 
ates from the expecte d values, consistent with previous analysis 
jLodato& Rice] [2004) suggesting that the local approximation is 
suitable when H/r < 0.1. 

4.2.1.3 The local and quasi-steady assumptions Figure \W\ 
also shows that, for gi n i t = 0.25, the temperature fluctuates by less 
than 10% and Q fluctuates by 10% - 20%, over the final 13 ORPs. 
This illustrates that all these discs settle into quasi-steady states that 
are marginally stable. The non-local transport fraction (Figure [TT) 



also remains low. The seemingly high £ for M* = O.5A/0 is due 
to its slightly elevated mass ratio in comparison to the other discs 
(Figure[8] bottom right panel). This, coupled with its comparatively 
lower sound speed and lower surface density (with the scale height 
kept constant) will boost the non-local transport fraction to a higher 
value than expected ab initio. However, its maximum value is still 
below that of the gi^it = 0.5 disc studied in this analysis (Simula- 
tion 2), so this is not inconsistent with expectations. 



4.2.2 The qi^it ~ 1 discs 

4.2.2.1 General Evolution Figure [T2l shows the profiles of the 
Qinit = 1 discs, averaged over the final 13 ORPs. The initial stellar 
masses are M* = 0.5M Q , M, = 1M Q , and M, = 2M . The 
discs grow hotter with increasing disc mass (with a flatter temper- 
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Figure 11. The non-local transport fraction for the qi n it = 0.25 simula- 
tions (Simulation 5 (solid line), Simulation 1 (dotted lines), Simulation 6 
(dashed lines) and Simulation 7 (dot-dashed lines)), averaged over the last 
13 ORPs. 

ature profile), while maintaining a similar surface density profile. 
This results in the higher disc mass simulations obtaining a flatter 
aspect ratio (top right panel). 

4.2.2.2 The a Approximation Figure[T3lshows that all the discs 
have similar qualitative a profiles, with a to tai being different to 
what would be expected if the local approximation were appropri- 
ate (a coo i). The value of the enhancement appears to increase with 
increasing disc mass, showing that while q dictates whether or not a 
disc deviates from the local approximation, the disc mass M<j con- 
trols the strength of this deviation (through its influence on E and 
ultimately the disc thickness). All three discs have aspect ratios in 
excess of 0. 1 for most of their radial extent, again consiste nt with 
previous predictions for non-locality jLodato & R ice 2004). 



local transport fraction dCossins et alj|200^) . moderate azimuthal 
Fourier mode amplitudes up to m ~ 8 (with increased power at 
m = 2), and maintain a str ictly self-regulated, quasi-steady state 
jLodato & Ricd 120041 l2005h . It has also been demonstrated that 
increasing stellar mass (while keeping q constant) does not sig- 
nificantly affect the efficacy of the viscous approximation, hold- 
ing over at least an order of magnitude in stellar mass. There is, 
however, some suggestion that there is some dependence on stel- 
lar mass with the M* = 5M© simulation showing some evidence 
for non-local transport corresponding to regions of the disc where 
H/r > 0.1. 

However, if the disc-to-star mass ratio qi n i t > 0.5, the az- 
imuthal m = 2 spiral modes begin to dominate. The strength of 
these global spiral waves introduces strong non-local torques, and 
are also subject to transient burst events. The disc stresses calcu- 
lated show that locally, in a time averaged sense, the amount of 
energy released through cooling does not match the thermal energy 
generated by the instability. It is likely that this excess energy is 
transported by the low-m mode global waves to larger (r > 40 au) 
radii where it can be lost through radiative cooling. This is a clear 
indication of globa l effects and is confir med by their high non-local 
transport fractions dCossins et al. 2009). Together, these violate the 
assumptions made to satisfy the viscous approximation. 

In summary, semi-analytic models are justified in using the 
viscous approximation to model realistic self-gravitating protostel- 
lar discs, provided that the parameter space studied does not in- 
clude discs that are too m assive, or geometrically thick. The c ur- 
rent semi-analytic models (Clarke 2009; Rice & Armitage 2009) in 
which the midplane cooling time is used to determine the effective 
gravitational a will, however, certainly underestimate the value of 
the effective viscosity in massive, geometrically thick discs. 



4.2.2.3 The local and quasi-steady assumptions Figure [14] 
also shows that the quasi-steady approximation also appears to be 
violated (Figure[T4l. The temperature fluctuates at values of ~ 20% 
and higher, with similar fluctuations in Q. The non-local transport 
fraction (bottom right panel in Figure [T3l > in all three cases is ~ 1 
or larger showing that the transport is very non-local. 
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5 CONCLUSIONS 

This work has studied in detail whether a local, viscous approxi- 
mation can accurately model the angular momentum transport in 
realistic, radiative, self-gravitating protostellar discs. For the vis- 
cous approximation to hold, the angul ar momentum transport must 
be local. If the analytical results of Gammiej J200 lh and others also 
hold (which calculate the stresses using the assumption that the dis- 
sipation rate matches the local cooling rate), the discs must also be 
in approximate thermodynamic equilibrium. 

A series of simulations using SPH with radiative transfer were 
carried out, and the effective viscosity generated by the gravita- 
tional instability was calculated directly from the Reynolds and 
gravitational stresses in the simulated discs. This was then com- 
pared with the expected viscosity, based on the assumption of local 
thermodynamic equilbrium, and the results analysed as a function 
of increasing disc-to-star mass ratio, and increasing stellar mass. 

The results show that if the discs have an initial disc-to-star 
mass ratio gi n i t < 0.5, and are geometrically thin (H/r ^ 0.1), the 
local viscous approximation performs well in calculating the angu- 
lar momentum transport. Such discs are shown to have a low non- 
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Figure 12. Azimuthally averaged radial profiles from the (jr init = 1 simulations (Simulation 8 (solid line), Simulation 3 (dotted lines), and Simulation 9 
(dashed lines)). The figures show the time average of each variable, taken from the last 13 ORPs. The top left panel shows the surface density profile, the 
top right shows the aspect ratio, the bottom left shows the midplane temperature, and the right hand panel shows the disc mass ratio q as a function of time. 
Artificial viscosity dominates inside 10 au, so data from inside this region is ignored. 
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Figure 13. The a parameter (Simulation 8 top left, Simulation 3 top right, & Simulation 9 bottom right), averaged over the last 13 ORPs of the simulations. 
The black line indicates the alpha calculated from Reynolds and gravitational stresses, the dashed line indicates the alpha calculated by the midplane cooling 
time at that radius, and the dotted line indicates the alpha calculated from the vertically averaged cooling time. The bottom right panel shows the non-local 
transport fraction for the <ji n i t = 1 simulations (Simulation 8 (solid line), Simulation 3 (dotted lines) and Simulation 9 (dashed lines)), averaged over the last 
13 ORPs. 
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Figure 14. Variation in the mean temperature profile (left) and the mean Toomre instability profile (right) for the q = 1 simulations (Simulation 8 (solid line), 
Simulation 3 (dotted lines), and Simulation 9 (dashed lines)), averaged over the last 13 ORPs. 
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